Abstract
Several regions show long-term changes in richness, generally increases. The rest of this study focuses on understanding these long-term changes through the lens of the geographic ranges of individual species and how geographic range relates to the colonization and extinction of those species, and thus to changes in richness.
***
rsType <- c("observed", "rarefied", "msom")[2]
rSMet_base <- c("observed"="propStrata", "rarefied"="range_size_samp", "msom"="range_size_mu")[rsType]
rSMet_annComm <- c("observed"="propStrata_avg", "rarefied"="range_size_samp_avg", "msom"="range_size_mu_avg")[rsType]
rSMet_ltComm <- c("observed"="propStrata_avg_ltAvg", "rarefied"="range_size_samp_avg_ltAvg", "msom"="range_size_mu_avg_ltAvg")[rsType]
rR <- c("observed"="naive_rich", "rarefied"=NA, "msom"="reg_rich")[3]
lR <- c("observed"="local_rich_obs", "rarefied"="local_rich_samp", "msom"="local_rich")[2]
betaD <- c("observed"="beta_div_obs", "msom"="beta_div_mu")[1]
There are four primary metrics:
- range size
- regional richness (gamma)
- local richness (alpha)
- spatial beta diversity
Each of these metrics can be based on raw observations or on MSOM estimates, and some can also be based on a rarefied (resampled, or ‘samp’) set of observations.
Furthermore, the range size metric can be based on a few combinations of cross-year and/or cross-species averages. The base metric for range size (no averaging) is used to to examine how range size changes with time until extinction and time after colonization. A cross-species average of a given year’s observed range sizes (_avg) is used for the rare species to see how this sub-group changes in its geographic distribution, or when making specialized comparisons. A cross-year average (ltAvg) is used to relate a species’ ‘typical’ range size to the total number of colonizations and extinctions it experienced during the time series. A cross-year (all years a species was observed) followed by a cross-species (species present in a given year) average (_avg_ltAvg) is used to compare species richness to the typical range size of species in the community.
# ============
# = Richness =
# ============
Number of species across all regions: 581
| reg | mu | sd | mu_sd |
|---|---|---|---|
| ai | 52.499 | 2.3143 | 5.299 |
| ebs | 103.611 | 5.5665 | 5.299 |
| gmex | 142.285 | 4.2663 | 5.299 |
| goa | 95.717 | 8.5061 | 5.299 |
| neus | 122.030 | 8.2371 | 5.299 |
| newf | 65.761 | 4.3653 | 5.299 |
| sa | 102.435 | 4.0355 | 5.299 |
| shelf | 46.300 | 3.2469 | 5.299 |
| wctri | 84.478 | 7.1532 | 5.299 |
# load regression results
load("../../pkgBuild/results/rich_naive_trend_kendall.RData")
rich_naive_trend_kendall[reg!="wcann",BH:=p.adjust(taup, method='BH')]
rich_naive_trend_kendall <- rich_naive_trend_kendall[
reg!="wcann",list(reg=reg, estimate=tau, BH=BH, pvalue=taup)
]
load("../../pkgBuild/results/rich_trend_kendall.RData")
rich_trend_kendall[reg!="wcann",BH:=p.adjust(pvalue, method="BH")]
rich_trend_kendall <- rich_trend_kendall[
reg!="wcann",list(reg=reg, estimate=tau, BH=BH, pvalue=pvalue)
]
# make plot
load("../results/rich_trend_kendall.RData") # load rich trend stats
richness_ts()
pdf("../../submitted_pdfs/final/Fig1.pdf", width=3.228, height=4.6114)
richness_ts()
dev.off()
quartz_off_screen 2
Figure 1. Time series of MSOM estimates of region richness. Each point is the posterior mean of regional richness in a year. Lines indicate long-term trends from fitted values of linear regression models predicting richness from time.
Looking for trends in species richness in both the naive estimates and the MSOM estimates. Using Kendall’s Tau_b, calculated for 1E4 resamplings of the posterior of richness. Tau is also calculated using a method that removes serial correlation to achieve independence of observations so that the p-values are correct.
In the Naive estimates, 7 regions had significant \(\tau_b\).
| reg | tau | Z | pvalue |
|---|---|---|---|
| ebs | 0.42141 | 3.31119 | 0.00093 |
| ai | 0.08918 | 0.37542 | 0.70735 |
| goa | 0.24549 | 1.08760 | 0.27677 |
| wctri | 0.60763 | 2.34639 | 0.01896 |
| wcann | 0.41618 | 1.67804 | 0.09334 |
| gmex | 0.09637 | 0.51845 | 0.60414 |
| sa | -0.21815 | -1.47174 | 0.14109 |
| neus | 0.19111 | 1.50292 | 0.13286 |
| shelf | 0.45244 | 4.06179 | 0.00005 |
| newf | 0.72656 | 3.80999 | 0.00014 |
In the MSOM estimates, 4 regions had significant \(\tau_{b}\).
Overall, ~half the regions show positive trends. No regions show significant negative trends (although SEUS is close, depending on the analysis). It’s a bit surprising how much removing the autocorrelation seemed to impact some of the trends (AI in particular, I think).
Manuscript paragraph:
Estimated slopes for long-term trends (Kendall’s \(\tau_b\)) in richness were positive for most regions. For Naive estimates, all the seven positive \(\tau\) were also significantly different from 0, whereas the two negative \(\tau_b\), Aleutian Islands and Southeast US, were not (Table S2). Long-term trends in MSOM estimates of richness had the same sign as Naive trends, except for Aleutian Islands, but now the only significant \(\tau_b\) were the following four regions: Eastern Bering Sea (\(\tau_b\) = 0.42141), West Coast US (\(\tau_b\) = 0.60763), Scotian Shelf (\(\tau_b\) = 0.45244 ), and Newfoundland (\(\tau_b\) = 0.72656 ) (Table 1). Richness trends were not significant in the Gulf of Mexico, Gulf of Alaska, and Northeast US, Aleutian Islands, Southeast US.
The number of species-regions: 863
The number of core species: 536
The number of both species: 263
The number of colonizing species: 60
The number of leaving species: 4
rangeSize_absenceTime(rSMet_base)
pdf("../../submitted_pdfs/final/Fig2.pdf", width=6.81102, height=3*(6.81102/6))
rangeSize_absenceTime(rSMet_base)
dev.off()
quartz_off_screen 2
print(rSMet_base)
rarefied
“range_size_samp”
Figure 4. Geographic range size vs years until extinction (A) and years after colonization (B). For visualization purposes, range size is averaged across species for each unique value on each axis, and a linear model fit through this average. Statistics in main text use unaggregated data. The horizontal axes were formulated as time until (since) the nearest upcoming (previous) absence. Because range size must be zero when either horizontal axis has a value of zero, points at (0,0) were excluded from figures and analyses.
Range size declines near an absence much more consistently than does range density; both are (relatively) low just before extinction and just after colonization. However, range density has much more variable intercepts among regions, whereas range size does not.
This makes sense, at least somewhat, because colonization and extinction events are defined at the site level; though the outcome isn’t necessitated by this formulation, because size could drop suddenly. In fact, when a species is absent, both its range size and its range density must be 0 (though, range density is technically calculated for only those sites that are occupied, so I supposed it’s technically undefined according to the equations I’m using).
I think the regressions for range size should omit an intercept, while the regressions for range density should have it. This might be hard to justify fully a priori (though see my thinking in previous paragraph), so I’ll probably just do a model selection and maybe discuess the difference if one model has an intercept and the other does not.
rangeTimeDT <- make_rangeTime(sizeName=rSMet_base)
# Set up -- region names and empty named lists
ur <- rangeTimeDT[,unique(reg)]
sTime_reg_mods0 <- structure(vector("list",length(ur)), .Names=ur)
sTime_reg_mods1 <- structure(vector("list",length(ur)), .Names=ur)
sTime_reg_mods2 <- structure(vector("list",length(ur)), .Names=ur)
sTime_reg_mods3 <- structure(vector("list",length(ur)), .Names=ur)
# Fit 3 models for each region
for(r in 1:length(ur)){
t_reg <- ur[r]
t_dat <- rangeTimeDT[reg==t_reg]
sTime_reg_mods0[[t_reg]] <- lme4::lmer(size ~ time + (time-1|spp), data=t_dat)
sTime_reg_mods1[[t_reg]] <- lme4::lmer(size ~ time + (time|spp), data=t_dat)
sTime_reg_mods2[[t_reg]] <- lme4::lmer(size ~ time + type + (time|spp), data=t_dat)
sTime_reg_mods3[[t_reg]] <- lme4::lmer(size ~ time * type + (time|spp), data=t_dat)
}
# Summarize each model
sTime_reg_smry0 <- smry_modList2(sTime_reg_mods0)
sTime_reg_smry1 <- smry_modList2(sTime_reg_mods1)
sTime_reg_smry2 <- smry_modList2(sTime_reg_mods2)
sTime_reg_smry3 <- smry_modList2(sTime_reg_mods3)
# Captions for each model
# Also a caption in case I want to report all models for all regions
capS0 <- sTime_reg_smry0[,unique(mod_call)]
capS1 <- sTime_reg_smry1[,unique(mod_call)]
capS2 <- sTime_reg_smry2[,unique(mod_call)]
capS3 <- sTime_reg_smry3[,unique(mod_call)]
timeRegMods_cap <- "Summary statistics for fits of predicting range size from years before extinction and years after colonization. These are mixed effect models of the form "
# Compare AIC values from each model fit
compAIC <- rbind(
cbind(sTime_reg_smry0[,list(reg,MargR2,CondR2,AIC)], mod=capS0),
cbind(sTime_reg_smry1[,list(reg,MargR2,CondR2,AIC)], mod=capS1),
cbind(sTime_reg_smry2[,list(reg,MargR2,CondR2,AIC)], mod=capS2),
cbind(sTime_reg_smry3[,list(reg,MargR2,CondR2,AIC)], mod=capS3)
)
setkey(compAIC, reg) # sort
# Get best model from each region, and from each region list best overall model
bestEach <- compAIC[,list("mod_call"=unique(mod)[which.min(AIC)]),by=c('reg')] # combinations of each region and its best model
bestOverall_name <- bestEach[,j={bm <- table(mod_call); names(bm)[which.max(bm)]}] # name of best overall model
bestOverall <- bestEach[,CJ(reg=reg,mod_call=bestOverall_name)] # get combos for each region w/ best overall model
bestEachOverall <- unique(rbind(bestEach,bestOverall)) # combinations for best overall and each
setkey(bestEachOverall, reg, mod_call) # sort
# get the worthy models
allSmry <- rbind(sTime_reg_smry0, sTime_reg_smry1, sTime_reg_smry2, sTime_reg_smry3, fill=TRUE) # combine results from all models and regions
bestModels <- allSmry[bestEachOverall, on=c('reg','mod_call')] # for each region, only select models that were best in that region or best overall
# remove columns that are NA for all rows b/c models w/ those parameters were never winners
loserNames <- sapply(bestModels, function(x)all(is.na(x))) # names of columns only pertaining to non-best models
loserNames <- names(loserNames)[loserNames] # okay, get the names for real, not just logic
bestModels[,c(loserNames):=list(NULL)] # drop the loser columns
| reg | mod_call | MargR2 | CondR2 | AIC | pval_time | BH_time | randGrp | Int | time | pval_type | BH_type | typepre_ext |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| ai | size~time+(time|spp) | 0.06079 | 0.6421 | -269.26 | 0.00576 | 0.006 | spp | 0.06754 | 0.00370 | NA | NA | NA |
| ebs | size~time+(time|spp) | 0.05742 | 0.8888 | -3766.58 | 0.00000 | 0.000 | spp | 0.04491 | 0.00447 | NA | NA | NA |
| gmex | size~time+(time-1|spp) | 0.08244 | 0.7549 | -1928.04 | 0.00076 | 0.001 | spp | 0.02208 | 0.00512 | NA | NA | NA |
| gmex | size~time+(time|spp) | NA | NA | -1924.04 | 0.00076 | 0.001 | spp | 0.02208 | 0.00512 | NA | NA | NA |
| goa | size~time+(time|spp) | 0.08819 | 0.6998 | -1394.63 | 0.00003 | 0.000 | spp | 0.02083 | 0.00222 | NA | NA | NA |
| neus | size~time+(time|spp) | 0.08595 | 0.7446 | -9567.84 | 0.00005 | 0.000 | spp | 0.01268 | 0.00141 | NA | NA | NA |
| newf | size~time+(time|spp) | 0.05885 | 0.8348 | -717.96 | 0.03647 | 0.036 | spp | 0.04496 | 0.00571 | NA | NA | NA |
| sa | size~time+(time|spp) | 0.12180 | 0.7792 | -2094.49 | 0.00215 | 0.003 | spp | 0.02733 | 0.00718 | NA | NA | NA |
| shelf | size~time+(time|spp) | 0.18780 | 0.6328 | -2646.81 | 0.00000 | 0.000 | spp | 0.03364 | 0.00247 | NA | NA | NA |
| shelf | size~time+type+(time|spp) | 0.19420 | 0.6391 | -2647.04 | 0.00000 | 0.000 | spp | 0.03762 | 0.00229 | 0.00058 | 0.005 | -0.01203 |
| wctri | size~time+(time|spp) | 0.14080 | 0.9525 | -471.45 | 0.00156 | 0.002 | spp | -0.00863 | 0.01116 | NA | NA | NA |
chosenModels <- bestModels[!is.na(MargR2) & mod_call!=("size~time+type+(time|spp)")]
mR2Range <- chosenModels[,range(MargR2)]
cR2Range <- chosenModels[,range(CondR2)]
pRange <- chosenModels[,range(pval_time)]
slopeRange <- chosenModels[,range(time)]
slopeAvg <- chosenModels[,mean(time)]
sd_ranef <- sapply(sTime_reg_mods1, function(x)sd(coef(x)[[1]][,'time']))
mean_ranef <- sapply(sTime_reg_mods1, function(x)mean(coef(x)[[1]][,'time']))
dRange_ranef <- sapply(sTime_reg_mods1, function(x)diff(range(coef(x)[[1]][,'time'])))
ranef_smry <- data.table(reg=names(sd_ranef), sd_ranef=sd_ranef, mean_ranef=mean_ranef, dRange_ranef=dRange_ranef)
ranef_smry[,c("sd_mean","dR_mean"):=list(sd_ranef/mean_ranef, dRange_ranef/mean_ranef)]
after correcting for multiple tests, all \(p \leq 0.03647; 0.05742 \leq mR^2 \leq 0.1878; 0.6328 \leq cR^2 \leq 0.9525\). For slope, \(1.40595 \leq \beta \leq 11.15633, \bar{\beta} = 4.82732\) percent per decade.
“Indeed, variation among species’ slopes was similar to the average slope(\(\hat{\sigma}_{\beta} = 6.47631\) percent per decade)” ###Conclusion for Range Change before Extinction/ after Colonization As stated with the larger (pooled regional data) regressions, range size responds more consistently/ strongly to an approaching extinction/ departure from colonization than does range density. Range size shrinks as extinction approaches, increases as time since colonization increases. In both cases, there was only 1 region for which the direction (into extinction, out of colonization) mattered: for range density, it was the Southeast US (effect of pre-extinction direction = 0.037, p=0.043 [BH correction], Table 12), and for range size it was Scotial Shelf (effect = -0.021, p = 0.002). Time was a significant predictor of range size in all regions; time was not a significant predictor of range density in Newfoundland (p = 0.846 [BH]), and the Scotial Shelf (p = 0.916).
These results support the hypothesis that species rarity is closely associated with proximity to extinction/ colonization, and therefore the probability that the species will contribute to a change in species richness. Furthermore, these relationships have not been directly quantified for entire assemblages. Competing hypotheses exist regarding how range should change approaching extinction and after colonization. We find that the change in range is similar regardless of direction. However, spatial scale was important. Although slopes were similar, the intercept for density was far larger than for range size — the fraction of tows containing a species in occupied sites may remain high even when extinction is imminent, and may similar be large even if colonization was recent. Range size, on the other hand, was much closer to 0 near an absence.
Overall, the results suggest that the spatial footprint of individual species is important for understanding changes in species richness. Furthermore, because species contributing most to the dynamics of richness are those that repeatedly colonize and go extinct, it is meaningful to look at a species’ long-term rarity in order to gauge whether it is likely to contribute to those long-term richness changes. Determining what drives the geographic range of individual species is probably a powerful way to anticipate richness changes.
transRange <- spp_master[,list(reg=reg, spp=spp, transient=(ce_categ!='neither'),rangeSize=range_size_samp)]
ur <- names(pretty_reg)
nr <- length(ur)
tR_mod <- structure(vector('list', length(ur)), .Names=ur)
for(r in 1:nr){
tR_mod[[ur[r]]] <- lmer(rangeSize~transient + (1|spp), data=transRange[reg==ur[r]])
}
tR_mod_smry <- smry_modList2(tR_mod)
kable(tR_mod_smry)
| reg | mod_call | MargR2 | CondR2 | AIC | pval_transient | BH_transient | randGrp | Int | transientTRUE |
|---|---|---|---|---|---|---|---|---|---|
| ai | rangeSize~transient+(1|spp) | 0.07447 | 0.8741 | -1233.9 | 0.02037 | 0.020 | spp | 0.25240 | -0.13085 |
| ebs | rangeSize~transient+(1|spp) | 0.16380 | 0.9404 | -8293.1 | 0.00000 | 0.000 | spp | 0.27140 | -0.19520 |
| gmex | rangeSize~transient+(1|spp) | 0.11860 | 0.8461 | -5264.1 | 0.00000 | 0.000 | spp | 0.19516 | -0.14126 |
| goa | rangeSize~transient+(1|spp) | 0.10370 | 0.9322 | -3308.5 | 0.00030 | 0.000 | spp | 0.17728 | -0.12959 |
| neus | rangeSize~transient+(1|spp) | 0.41200 | 0.8827 | -11018.4 | 0.00000 | 0.000 | spp | 0.19523 | -0.17705 |
| newf | rangeSize~transient+(1|spp) | 0.08736 | 0.9014 | -2573.9 | 0.00214 | 0.002 | spp | 0.18857 | -0.12827 |
| sa | rangeSize~transient+(1|spp) | 0.28070 | 0.8110 | -3536.0 | 0.00000 | 0.000 | spp | 0.34518 | -0.28582 |
| shelf | rangeSize~transient+(1|spp) | 0.30740 | 0.8681 | -4145.8 | 0.00000 | 0.000 | spp | 0.28011 | -0.22187 |
| wctri | rangeSize~transient+(1|spp) | 0.11340 | 0.8470 | -1328.3 | 0.00004 | 0.000 | spp | 0.24779 | -0.17944 |
transDiffMean <- tR_mod_smry[,mean(transientTRUE)]*100
sign_tDM <- sign(transDiffMean)
In mixed effects models with intecepts varying among species and transient versus core as a categorical predictor, transient species had range sizes that were 17.65964 (% occupancy) smaller than the ranges of core species (average intercept adjustment; all \(p \leq\) 0.02 after correcting for multiple tests).
rich_geoRange(rSMet_ltComm, leg=TRUE, legPan=1, panLab=FALSE)
pdf("../../submitted_pdfs/final/Fig3.pdf", width=3.228, height=3.228)
rich_geoRange(rSMet_ltComm, leg=TRUE, legPan=1, panLab=FALSE)
dev.off()
quartz_off_screen 2
Figure 3. Species richness vs geographic range size. Range size is presented as each species’ long-term average of the proportion of sites it occupied. Solid lines are linear regressions with MSOM richness as the response and the horizontal axis and an intercept as the predictors.
Range size is a pretty good predictor of species richness. I think I had originally missed the range size relationship b/c I hadn’t done the same aggregating procedure. The interpretation I have is that richness is highest when you have a bunch of rare species.
range_reg <- make_range_reg(dens="propTow_occ_avg", size=rSMet_ltComm)
# Fit different models to the whole data set
# Models vary in which/ how parameters vary among regions
rSize_mods <- list()
rSize_mods[[1]] <- lm(rich ~ size, data=range_reg) # simple
rSize_mods[[2]] <- lm(rich ~ size*reg, data=range_reg) # slope and intercept ind. among regs
rSize_mods[[3]] <- lme4::lmer(rich ~ size + (1|reg), data=range_reg) # intercept varies randomly among regs
rSize_mods[[4]] <- lme4::lmer(rich ~ size + (size|reg), data=range_reg) # slope & int. vary randomly among regs
rich_size_smry <- smry_modList(rSize_mods) # summarize to get R^2, AIC, p-vals, etc
## Note: method with signature 'sparseMatrix#ANY' chosen for function 'kronecker',
## target signature 'dgCMatrix#ngCMatrix'.
## "ANY#sparseMatrix" would also be valid
# Fit simple model to each region separately
ur <- range_reg[,unique(reg)]
rSize_reg_mods <- structure(vector("list",length(ur)), .Names=ur)
for(r in 1:length(ur)){
rSize_reg_mods[[r]] <- lm(rich ~ size, data=range_reg[reg==ur[r]])
}
rich_s_reg_smry <- data.table(reg=ur, smry_modList(rSize_reg_mods)) # summarize
setnames(rich_s_reg_smry, old=c("Marginal","Conditional","p.value"), new=c("MargR2","CondR2","pval")) # rename
setkey(rich_s_reg_smry, reg, Class, mod_call, predictor) # sort
rich_s_reg_smry[, BH:=round(p.adjust(pval, "BH"), 3), by=c("mod_call", "predictor")] # correct p-vals for multiple tests
smry_formula <- formula(reg+Class+mod_call+MargR2+CondR2+AIC~predictor) # formula to organize
rich_s_reg_smry <- dcast(rich_s_reg_smry, smry_formula, value.var=c("pval", "BH")) # organize; helpful if multiple predictors
| reg | Class | mod_call | MargR2 | CondR2 | AIC | pval_size | BH_size |
|---|---|---|---|---|---|---|---|
| ai | lm | rich ~ size | 0.7348 | NA | 43.222 | 0.0004 | 0.000 |
| ebs | lm | rich ~ size | 0.8696 | NA | 136.247 | 0.0000 | 0.000 |
| gmex | lm | rich ~ size | 0.7684 | NA | 77.674 | 0.0000 | 0.000 |
| goa | lm | rich ~ size | 0.7908 | NA | 77.173 | 0.0000 | 0.000 |
| neus | lm | rich ~ size | 0.8609 | NA | 167.633 | 0.0000 | 0.000 |
| newf | lm | rich ~ size | 0.9151 | NA | 58.069 | 0.0000 | 0.000 |
| sa | lm | rich ~ size | 0.3582 | NA | 134.595 | 0.0016 | 0.002 |
| shelf | lm | rich ~ size | 0.8999 | NA | 123.548 | 0.0000 | 0.000 |
| wctri | lm | rich ~ size | 0.9510 | NA | 42.516 | 0.0000 | 0.000 |
| mod_call | reg | Class | MargR2 | CondR2 | AIC | pval_size | BH_size |
|---|---|---|---|---|---|---|---|
| rich ~ size | NA | NA | 0.7943 | NA | 95.631 | 0.00022 | 0.00022 |
Manuscript Text:
Furthermore, species richness was negatively correlated with the community range index (CRI) in each region (Fig. 3; separate linear regression for each region, for slope all corrected \(p \leq\) 0.002, \(0.3582 \leq R^2 \leq 0.951\), average \(R^2 = 0.7943\) ).
**Other Notes:**
All models are pretty good predictors. Well, the most basic model kinda sucks I guess. It needs to account for some of the between-region variation.
Predicting Richness: Range Size or Density? As far as picking one or the other, it doesn’t end up mattering much. Range size is a lot better than density in NEUS, and density outperforms size in AI. Otherwise, size as a slight edge over density on average, although both predictors are significant in all regions.
Conclusion for Richness and Geographic Range The two metrics of geographic range are well correlated. Furthermore, richness can be predicted pretty well using regressions with either as a predictor. There are large differences among regions, though. This is probably because richness is not readily comparable among most regions. Regions vary mostly in their intercept values, and they have fairly similar slopes (though they are not identical, and model fits improve when allowing slopes to vary among regions; it’s just that the improvement is small compared to allowing intercepts to vary among regions).
The interpretation of the result that geographic distribution predicts species richness is likely associated with species rarity. When the average range density or range size of a community is low, it means it has a lot of species that are rare (at either spatial scale). It’s these rare species that come and go, and form the dynamics of richness that we observe. When that dynamical value is high, it implies that an above-average number of the dynamic species are present. Because those species are transient (dynamic), they are also probably rare.
eval(neitherPres_subColor())
plot_rangeSize_FullTrans(range_type="range_size_samp")
pdf("../../submitted_pdfs/final/Fig4.pdf", width=6.81102, height=6.81102)
plot_rangeSize_FullTrans(range_type="range_size_samp")
dev.off()
quartz_off_screen 2
transMod <- structure(vector("list",length(ur)), .Names=ur)
coreMod <- structure(vector("list",length(ur)), .Names=ur)
ctMod <- structure(vector("list",length(ur)), .Names=ur)
ctMod2 <- ctMod
for(r in 1:length(ur)){
rr <- ur[r]
td <- spp_master[eval(rE_logic1)]
td[,year:=year-min(year)]
transMod[[r]] <- lmer(range_size_samp~year+(1|spp), data=td)
td2 <- spp_master[eval(rE_logic3)]
td2[,year:=year-min(year)]
coreMod[[r]] <- lmer(range_size_samp~year+(1|spp), data=td2)
td3 <- spp_master[reg==rr & present==1]
td3[,year:=year-min(year)]
td3[,categ2:=c("core","transient")[(ce_categ!="neither")+1L]]
ctMod[[r]] <- lmer(range_size_samp~year*categ2+(1|spp), data=td3)
ctMod2[[r]] <- lm(range_size_samp~year*categ2, data=td3)
# ctMod2[[r]] <- lmer(range_size_samp~year+(year|spp), data=td)
# ctMod2[[r]] <- lmer(range_size_samp~year+(year|spp), data=td3)
# ctMod2[[r]] <- lmer(range_size_samp~year*categ2+(year|spp), data=td3)
# td3[,auto.arima(range_size_samp, xreg=as.matrix(cbind(year=year, categ2=(categ2=="transient"), year_categ2=(year*(categ2=="transient")))), allowdrift=FALSE, d=0)]
}
transExpansionStats <- smry_modList2(transMod)
coreExpansionStats <- smry_modList2(coreMod)
ctExpansionStats <- smry_modList2(ctMod)
# ctExpansionStats2 <- smry_modList2(ctMod2)
meanNum <- function(x, ...){
if(storage.mode(x) %in% c("double", "integer", "logical")){
mean(x, ...)
}else{
NA
}
}
kable(rbind(transExpansionStats,lapply(transExpansionStats, meanNum), transExpansionStats[BH_year<0.05, lapply(.SD, meanNum)][,reg:="signifYear"]))
| reg | mod_call | MargR2 | CondR2 | AIC | pval_year | BH_year | randGrp | Int | year |
|---|---|---|---|---|---|---|---|---|---|
| ai | range_size_samp~year+(1|spp) | 0.04978 | 0.63660 | -266.83 | 1e-05 | 0 | spp | 0.05853 | 0.00343 |
| ebs | range_size_samp~year+(1|spp) | 0.01409 | 0.82480 | -3383.90 | 0e+00 | 0 | spp | 0.05250 | 0.00143 |
| gmex | range_size_samp~year+(1|spp) | 0.03786 | 0.63360 | -1543.33 | 0e+00 | 0 | spp | 0.02259 | 0.00363 |
| goa | range_size_samp~year+(1|spp) | 0.07234 | 0.57180 | -1322.23 | 0e+00 | 0 | spp | 0.01342 | 0.00186 |
| neus | range_size_samp~year+(1|spp) | 0.02210 | 0.54520 | -8917.73 | 0e+00 | 0 | spp | 0.01091 | 0.00044 |
| newf | range_size_samp~year+(1|spp) | 0.03951 | 0.66150 | -719.18 | 0e+00 | 0 | spp | 0.03018 | 0.00315 |
| sa | range_size_samp~year+(1|spp) | 0.05840 | 0.34100 | -1760.63 | 0e+00 | 0 | spp | 0.08827 | -0.00252 |
| shelf | range_size_samp~year+(1|spp) | 0.12960 | 0.56570 | -2508.97 | 0e+00 | 0 | spp | 0.01941 | 0.00170 |
| wctri | range_size_samp~year+(1|spp) | 0.11650 | 0.45570 | -332.63 | 0e+00 | 0 | spp | -0.01374 | 0.00488 |
| NA | NA | 0.06002 | 0.58177 | -2306.16 | 0e+00 | 0 | NA | 0.03134 | 0.00200 |
| signifYear | NA | 0.06002 | 0.58177 | -2306.16 | 0e+00 | 0 | NA | 0.03134 | 0.00200 |
kable(rbind(coreExpansionStats,lapply(coreExpansionStats, meanNum), coreExpansionStats[BH_year<0.05, lapply(.SD, meanNum)][,reg:="signifYear"]))
| reg | mod_call | MargR2 | CondR2 | AIC | pval_year | BH_year | randGrp | Int | year |
|---|---|---|---|---|---|---|---|---|---|
| ai | range_size_samp~year+(1|spp) | 0.00004 | 0.90040 | -969.63 | 0.65562 | 0.65600 | spp | 0.24991 | 0.00015 |
| ebs | range_size_samp~year+(1|spp) | 0.00289 | 0.94840 | -5170.31 | 0.00000 | 0.00000 | spp | 0.24800 | 0.00156 |
| gmex | range_size_samp~year+(1|spp) | 0.00047 | 0.84690 | -3827.05 | 0.01919 | 0.02500 | spp | 0.18865 | 0.00081 |
| goa | range_size_samp~year+(1|spp) | 0.00084 | 0.94090 | -2142.47 | 0.00063 | 0.00100 | spp | 0.16622 | 0.00070 |
| neus | range_size_samp~year+(1|spp) | 0.01428 | 0.85840 | -4630.15 | 0.00000 | 0.00000 | spp | 0.16285 | 0.00209 |
| newf | range_size_samp~year+(1|spp) | 0.00150 | 0.91320 | -1913.24 | 0.00024 | 0.00000 | spp | 0.20137 | -0.00171 |
| sa | range_size_samp~year+(1|spp) | 0.00268 | 0.77120 | -2162.97 | 0.00001 | 0.00000 | spp | 0.36634 | -0.00176 |
| shelf | range_size_samp~year+(1|spp) | 0.00041 | 0.83910 | -2177.36 | 0.09164 | 0.10300 | spp | 0.28729 | -0.00036 |
| wctri | range_size_samp~year+(1|spp) | 0.02543 | 0.89720 | -1165.98 | 0.00000 | 0.00000 | spp | 0.18834 | 0.00440 |
| NA | NA | 0.00540 | 0.87952 | -2684.35 | 0.08526 | 0.08722 | NA | 0.22878 | 0.00065 |
| signifYear | NA | 0.00687 | 0.88231 | -3001.74 | 0.00287 | 0.00371 | NA | 0.21740 | 0.00087 |
kable(rbind(ctExpansionStats,lapply(ctExpansionStats, meanNum), ctExpansionStats[BH_year<0.05, lapply(.SD, meanNum)][,reg:="signifYear"], ctExpansionStats[`BH_year:categ2`<0.05, lapply(.SD, meanNum)][,reg:="signifYear:categ2"]))
| reg | mod_call | MargR2 | CondR2 | AIC | pval_categ2 | pval_year | pval_year:categ2 | BH_categ2 | BH_year | BH_year:categ2 | randGrp | Int | year | categ2transient | year:categ2transient |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| ai | range_size_samp~year*categ2+(1|spp) | 0.07939 | 0.87930 | -1227.8 | 0.01933 | 0.00886 | 0.00002 | 0.01900 | 0.01000 | 0.00000 | spp | 0.24991 | 0.00015 | -0.19203 | 0.00330 |
| ebs | range_size_samp~year*categ2+(1|spp) | 0.16860 | 0.94380 | -8437.4 | 0.00000 | 0.00000 | 0.56113 | 0.00000 | 0.00000 | 0.59200 | spp | 0.24800 | 0.00156 | -0.19537 | -0.00014 |
| gmex | range_size_samp~year*categ2+(1|spp) | 0.12140 | 0.84870 | -5271.7 | 0.00000 | 0.00000 | 0.00008 | 0.00000 | 0.00000 | 0.00000 | spp | 0.18865 | 0.00081 | -0.16610 | 0.00282 |
| goa | range_size_samp~year*categ2+(1|spp) | 0.10660 | 0.93470 | -3318.8 | 0.00023 | 0.00000 | 0.00259 | 0.00000 | 0.00000 | 0.00400 | spp | 0.16622 | 0.00070 | -0.15278 | 0.00116 |
| neus | range_size_samp~year*categ2+(1|spp) | 0.42250 | 0.89330 | -11307.9 | 0.00000 | 0.00000 | 0.00000 | 0.00000 | 0.00000 | 0.00000 | spp | 0.16285 | 0.00209 | -0.15195 | -0.00165 |
| newf | range_size_samp~year*categ2+(1|spp) | 0.09326 | 0.90410 | -2570.9 | 0.00239 | 0.04512 | 0.00000 | 0.00300 | 0.04500 | 0.00000 | spp | 0.20137 | -0.00171 | -0.17034 | 0.00481 |
| sa | range_size_samp~year*categ2+(1|spp) | 0.28460 | 0.81490 | -3550.6 | 0.00000 | 0.00000 | 0.22953 | 0.00000 | 0.00000 | 0.29500 | spp | 0.36634 | -0.00176 | -0.27804 | -0.00079 |
| shelf | range_size_samp~year*categ2+(1|spp) | 0.31270 | 0.87250 | -4165.8 | 0.00000 | 0.00903 | 0.00000 | 0.00000 | 0.01000 | 0.00000 | spp | 0.28729 | -0.00036 | -0.26814 | 0.00207 |
| wctri | range_size_samp~year*categ2+(1|spp) | 0.14520 | 0.87880 | -1468.3 | 0.00001 | 0.00000 | 0.59204 | 0.00000 | 0.00000 | 0.59200 | spp | 0.18834 | 0.00440 | -0.20238 | 0.00047 |
| NA | NA | 0.19269 | 0.88557 | -4591.0 | 0.00244 | 0.00700 | 0.15393 | 0.00244 | 0.00722 | 0.16478 | NA | 0.22878 | 0.00065 | -0.19746 | 0.00134 |
| signifYear | NA | 0.19269 | 0.88557 | -4591.0 | 0.00244 | 0.00700 | 0.15393 | 0.00244 | 0.00722 | 0.16478 | NA | 0.22878 | 0.00065 | -0.19746 | 0.00134 |
| signifYear:categ2 | NA | 0.18931 | 0.88877 | -4643.8 | 0.00366 | 0.01050 | 0.00045 | 0.00367 | 0.01083 | 0.00067 | NA | 0.20938 | 0.00028 | -0.18356 | 0.00209 |
Examining how range size changed over time for “core” and “transient” groups of species. Core species are those that are present in all years. Transient species are those that are not present in all years. 9 regions had significant trends in range size for transient species (average slope for regions with significant trend = 0.002), and 7 regions had significant trends in range size for core species (ai, shelf did not have significant trends for core species; averge slope for regions with significant trend = 8.7071910^{-4}).
Given that many regions showed positive trends in range size for both species groups, it is prudent to compare trends between these groups, to see if they differ. To this end, I analyzed the groups together, and tested the significance of an interaction term that adjusted the a baseline slope if the species was in the transient group.
Note that transient species are very much expected to have a smaller range size than the core group, particularly at the beginning of the time series; I’m not dynamically reporting the significance of the difference in intercepts here, but at the time of writing all regions showed significant differences in intercept between groups, with the transient group having a smaller intercept than the core group.
Interestingly, 9 regions had significant trends (main effect), and newf, sa, shelf had significant negative trends, while ai, ebs, gmex, goa, neus, wctri had significant positive trends. ebs, sa, wctri were the regions without significant interactions – (at the time of this writing) these are not the same regions that lacked significant trends in the core group. If we sum the main effect and transient-interaction together (ignoring significance), we see that sa was/were the only region/s with negative estimated trends for transient species; all others had positive trends in range for transients. Some regions had significant trends for both core and transients, but the direction of these trends differed between the two groups: newf, shelf had core species with trends in range size that were significantly less than 0 but transient species with trends in range size that were significantly greater than 0.
Southeast U.S. core species long-term slope \(\beta_Y=\)-1.76278 (percent per decade)
Newfoundland core species slope \(\beta_Y=\)-1.7062 (percent per decade)
Scotian Shelf core species slope \(\beta_Y=\)-0.35899 (percent per decade)
For core species with positive long-term slopes, the average was **${}_{Y+}=\(1.61859** (percent per decade). After correcting for multipel testing, all **\)p $ 0.045**
However, the interaction term indicated that the slopes of core and transient species were different in 6 of the regions (\(p \leq\) 0.004), including positive interactions in Newfoundland (\(\beta_{Y \times T} =\) 4.81063) and Scotian Shelf (\(\beta_{Y \times T} =\) 2.0651). Excluding the Southeast U.S., the average long-term trend in range size for transient species was 2.56144
plotSiteMap()
spp_master[,j={cat('\n\n\n\n\n');cat('**',reg[1],'**',sep=''); cat('\n\n'); cat(paste(sort(unique(spp)),collapse=', '));cat('\n');NULL},by='reg']
ai
Anoplopoma fimbria, Aptocyclus ventricosus, Atheresthes evermanni, Atheresthes stomias, Bathymaster signatus, Bathyraja aleutica, Bathyraja maculata, Bathyraja parmifera, Bathyraja taranetzi, Berryteuthis magister, Careproctus rastrinus, Cheiraster dawsoni, Chionoecetes bairdi, Chlamys rubida, Dasycottus setiger, Enteroctopus dofleini, Erimacrus isenbeckii, Eumicrotremus orbis, Fusitriton oregonensis, Gadus macrocephalus, Glyptocephalus zachirus, Gorgonocephalus eucnemis, Gymnocanthus galeatus, Hemilepidotus jordani, Hemilepidotus zapus, Hemitripterus bolini, Hippoglossoides elassodon, Hippoglossus stenolepis, Hyas lyratus, Lethotremus muticus, Lithodes aequispinus, Malacocottus zonurus, Microstomus pacificus, Modiolus modiolus, Myoxocephalus polyacanthocephalus, Oregonia gracilis, Pagurus aleuticus, Pandalus eous, Pandalus tridens, Pleurogrammus monopterygius, Podothecus accipenserinus, Reinhardtius hippoglossoides, Rossia pacifica, Sarritor frenatus, Sebastes alutus, Sebastes borealis, Sebastes polyspinis, Sebastolobus alascanus, Strongylocentrotus droebachiensis, Theragra chalcogramma, Triglops forficatus, Triglops macellus, Triglops scepticus, Zaprora silenus
ebs
Aforia circinata, Ammodytes hexapterus, Argis dentata, Argis lar, Artediellus pacificus, Asterias amurensis, Atheresthes stomias, Bathymaster signatus, Beringius behringi, Boltenia ovifera, Buccinum angulosum, Buccinum plectrum, Buccinum polare, Buccinum scalariforme, Chionoecetes bairdi, Chionoecetes opilio, Chrysaora melanaster, Ciliatocardium ciliatum, Clinopegma magnum, Clupea pallasii pallasii, Crangon dalli, Crossaster papposus, Cucumaria fallax, Cyanea capillata, Dasycottus setiger, Echinarachnius parma, Eleginus gracilis, Erimacrus isenbeckii, Eunoe depressa, Eunoe nodosa, Euspira pallida, Evasterias echinosoma, Fusitriton oregonensis, Gadus macrocephalus, Gersemia rubiformis, Glebocarcinus oregonensis, Glyptocephalus zachirus, Gorgonocephalus eucnemis, Grandicrepidula grandis, Halichondria panicea, Halocynthia aurantium, Hemilepidotus jordani, Hemitripterus bolini, Hexagrammos stelleri, Hiatella arctica, Hippoglossoides elassodon, Hippoglossoides robustus, Hippoglossus stenolepis, Hyas coarctatus, Hyas lyratus, Icelus spatula, Icelus spiniger, Leptasterias arctica, Leptasterias groenlandica, Leptasterias polaris, Leptoclinus maculatus, Lethasterias nanimensis, Limanda aspera, Limanda proboscidea, Limanda sakhalinensis, Liparis gibbus, Liponema brevicorne, Lycodes brevipes, Lycodes palearis, Lycodes raridens, Mactromeris polynyma, Mallotus villosus, Metridium farcimen, Musculus discors, Neocrangon communis, Neptunea borealis, Neptunea heros, Neptunea lyrata, Neptunea pribiloffensis, Neptunea ventricosa, Occella dodecaedron, Ophiura sarsii, Oregonia gracilis, Pandalus eous, Pandalus goniurus, Paralithodes camtschaticus, Paralithodes platypus, Patinopecten caurinus, Platichthys stellatus, Pleuronectes quadrituberculatus, Plicifusus kroeyeri, Podothecus accipenserinus, Pteraster obscurus, Pyrulofusus deformis, Pyrulofusus melonis, Reinhardtius hippoglossoides, Rhamphostomella costata, Sarritor frenatus, Serratiflustra serrulata, Serripes groenlandicus, Siliqua alta, Stomphia coccinea, Strongylocentrotus droebachiensis, Styela rustica, Tellina lutea, Telmessus cheiragonus, Thaleichthys pacificus, Theragra chalcogramma, Trichodon trichodon, Triglops pingelii, Triglops scepticus, Tritonia diomedea, Urticina crassicornis, Volutopsius fragilis, Volutopsius middendorffii
gmex
Achelous spinicarpus, Achelous spinimanus, Amusium papyraceum, Anadara baughmani, Anasimus latus, Anchoa hepsetus, Anchoa lyolepis, Anchoa mitchilli, Anchoviella perfasciata, Ancylopsetta dilecta, Ancylopsetta ommata, Arenaeus cribrarius, Ariopsis felis, Balistes capriscus, Bellator militaris, Bollmannia communis, Bregmaceros atlanticus, Brevoortia patronus, Brotula barbata, Calappa sulcata, Callinectes sapidus, Callinectes similis, Caranx crysos, Caulolatilus intermedius, Centropristis philadelphica, Chaetodipterus faber, Chloroscombrus chrysurus, Chrysaora quinquecirrha, Citharichthys macrops, Citharichthys spilopterus, Cyclopsetta chittendeni, Cynoscion arenarius, Cynoscion nothus, Decapterus punctatus, Diplectrum bivittatum, Diplectrum formosum, Doryteuthis pleii, Engraulis eurystole, Engyophrys senta, Etropus crossotus, Etropus microstomus, Etrumeus teres(sadina), Eucinostomus gula, Farfantepenaeus aztecus, Farfantepenaeus duorarum, Fowlerichthys radiosus, Gymnachirus texae, Gymothorax igromargiatus, Haemulon aurolineatum, Halieutichthys aculeatus, Harengula jaguana, Hepatus epheliticus, Hoplunnis macrurus, Kathetostoma albigutta, Lagocephalus laevigatus, Lagodon rhomboides, Larimus fasciatus, Leiostomus xanthurus, Lepophidium brevibarbe, Lepophidium jeannae, Libinia dubia, Libinia emarginata, Litopenaeus setiferus, Loligo pealeii, Lolliguncula brevis, Lutjanus campechanus, Lutjanus synagris, Menticirrhus americanus, Micropogonias undulatus, Mullus auratus, Mustelus canis, Neobythites gilli, Neomerinthe hemingwayi, Ophidion holbrookii, Ophidion welshi, Opisthonema oglinum, Orthopristis chrysoptera, Ovalipes floridaus, Ovalipes stephensoni, Paralichthys lethostigma, Parapenaeus politus, Pareques umbrosus, Parthenopoides massena, Peprilus burti, Peprilus paru, Pericoptus punctatus, Peristedion gracile, Polydactylus octoemus, Pontinus longispinis, Porichthys plectrodon, Portunus gibbesii, Priacanthus arenatus, Prionotus longispinosus, Prionotus ophryas, Prionotus roseus, Prionotus rubio, Prionotus scitulus, Prionotus stearnsi, Prionotus tribulus, Pristipomoides aquilonaris, Raioides louisiaesis, Raja texana, Rhizoprionodon terraenovae, Rhomboplites aurorubens, Rhynchoconger flavus, Rimapenaeus constrictus, Rimapenaeus similis, Sardinella aurita, Saurida brasiliensis, Saurida caribbaea, Scomber japonicus, Scomberomorus cavalla, Scomberomorus maculatus, Scorpaena calcarata, Selar crumenophthalmus, Selene setapinnis, Seriola dumerili, Serranus atrobranchus, Sicyonia brevirostris, Sicyonia burkenroadi, Sicyonia dorsalis, Soleocera vioscai, Sphoeroides parvus, Sphyraena guachancho, Squatina dumeril, Squilla chydaea, Squilla empusa, Steindachneria argentea, Stellifer lanceolatus, Stenotomus caprinus, Stephanolepis hispidus, Syacium gunteri, Syacium papillosum, Symphurus civitatum, Symphurus diomedeanus, Symphurus plagiusa, Synodus foetens, Synodus poeyi, Trachurus lathami, Trichiurus lepturus, Trichopsetta ventralis, Upeneus parvus, Urophycis cirrata, Urophycis floridana
goa
Albatrossia pectoralis, Anoplopoma fimbria, Aphrocallistes vastus, Aptocyclus ventricosus, Arctomelon stearnsii, Atheresthes stomias, Bathymaster signatus, Bathyraja aleutica, Bathyraja interrupta, Bathyraja parmifera, Berryteuthis magister, Brisaster latifrons, Chionoecetes bairdi, Chlamys rubida, Clupea pallasii pallasii, Ctenodiscus crispatus, Cucumaria fallax, Cyanea capillata, Dasycottus setiger, Doris odhneri, Enteroctopus dofleini, Eumicrotremus orbis, Eumicrotremus phrynoides, Fusitriton oregonensis, Gadus macrocephalus, Glebocarcinus oregonensis, Glyptocephalus zachirus, Gorgonocephalus eucnemis, Gymnocanthus galeatus, Halocynthia aurantium, Hemilepidotus jordani, Hemitripterus bolini, Hexagrammos decagrammus, Hippasteria phrygiana, Hippoglossoides elassodon, Hippoglossus stenolepis, Isopsetta isolepis, Labidochirus splendescens, Leuroglossus schmidti, Limanda aspera, Lopholithodes foraminatus, Lumpenella longirostris, Lycodes brevipes, Lycodes diapterus, Lycodes palearis, Lyopsetta exilis, Malacocottus zonurus, Mallotus villosus, Mesocentrotus franciscanus, Microgadus proximus, Microstomus pacificus, Modiolus modiolus, Myoxocephalus jaok, Myoxocephalus polyacanthocephalus, Neptunea lyrata, Neptunea pribiloffensis, Oncorhynchus keta, Oncorhynchus tshawytscha, Ophiodon elongatus, Oregonia gracilis, Pagurus aleuticus, Pagurus brandti, Pagurus ochotensis, Pandalopsis dispar, Pandalus eous, Pandalus jordani, Parophrys vetulus, Pasiphaea pacifica, Patinopecten caurinus, Platichthys stellatus, Pleurogrammus monopterygius, Pleuronectes quadrituberculatus, Podothecus accipenserinus, Psettichthys melanostictus, Pycnopodia helianthoides, Raja binoculata, Raja rhina, Rhinolithodes wosnessenskii, Sebastes alutus, Sebastes babcocki, Sebastes borealis, Sebastes polyspinis, Sebastes ruberrimus, Sebastes variegatus, Sebastes zacentrus, Sebastolobus alascanus, Somniosus pacificus, Squalus acanthias, Strongylocentrotus droebachiensis, Strongylocentrotus fragilis, Strongylocentrotus pallidus, Stylasterias forreri, Thaleichthys pacificus, Theragra chalcogramma, Trichodon trichodon, Triglops forficatus, Triglops scepticus, Zaprora silenus
neus
Acipenser oxyrinchus, Alosa aestivalis, Alosa pseudoharengus, Alosa sapidissima, Amblyraja radiata, Ammodytes dubius, Anarhichas lupus, Anchoa hepsetus, Anchoa mitchilli, Antigonia capros, Argentina silus, Argentina striata, Aspidophoroides monopterygius, Atlantopandalus propinqvus, Bairdiella chrysoura, Bathynectes longispina, Bathypolypus arcticus, Brevoortia tyrannus, Brosme brosme, Cancer borealis, Cancer plebejus, Carcharhinus plumbeus, Centropristis philadelphica, Centropristis striata, Chaceon quinquedens, Chionoecetes opilio, Chlorophthalmus agassizi, Citharichthys arctifrons, Clupea harengus, Conger oceanicus, Crangon septemspinosa, Cryptacanthodes maculatus, Cyclopterus lumpus, Cynoscion regalis, Dasyatis centroura, Dasyatis say, Decapterus punctatus, Dichelopandalus leptocerus, Dipturus laevis, Enchelyopus cimbrius, Etropus microstomus, Etrumeus teres(sadina), Foetorepus agassizii, Gadus morhua, Glyptocephalus cynoglossus, Gymnura altavela, Helicolenus dactylopterus, Hemitripterus americanus, Hippoglossina oblonga, Hippoglossoides platessoides, Hippoglossus hippoglossus, Homarus americanus, Illex illecebrosus, Lagodon rhomboides, Larimus fasciatus, Lebbeus polaris, Leiostomus xanthurus, Lepophidium profundorum, Leptoclinus maculatus, Leucoraja erinacea, Leucoraja garmani, Leucoraja ocellata, Limanda ferruginea, Limulus polyphemus, Liparis atlanticus, Lithodes maja, Loligo pealeii, Lolliguncula brevis, Lophius americanus, Lopholatilus chamaeleonticeps, Lumpenus lampretaeformis, Lycenchelys verrillii, Macrorhamphosus scolopax, Malacoraja senta, Maurolicus weitzmani, Melanogrammus aeglefinus, Melanostigma atlanticum, Menidia menidia, Menticirrhus americanus, Menticirrhus saxatilis, Merluccius albidus, Merluccius bilinearis, Micropogonias undulatus, Monolene sessilicauda, Morone saxatilis, Mustelus canis, Myliobatis freminvillii, Myoxocephalus aenaeus, Myoxocephalus octodecemspinosus, Myxine glutinosa, Nemichthys scolopaceus, Ophidion marginatum, Osmerus mordax, Ovalipes ocellatus, Ovalipes stephensoni, Pandalus borealis, Pandalus montagui, Paralichthys dentatus, Parasudis truculenta, Pasiphaea multidentata, Peprilus triacanthus, Peristedion miniatum, Petromyzon marinus, Phycis chesteri, Placopecten magellanicus, Pollachius virens, Polymixia lowei, Pomatomus saltatrix, Pontinus longispinis, Pontophilus norvegicus, Prionotus alatus, Prionotus carolinus, Prionotus evolans, Pseudopleuronectes americanus, Raja eglanteria, Scomber scombrus, Scophthalmus aquosus, Scyliorhinus retifer, Sebastes fasciatus, Sicyonia brevirostris, Squalus acanthias, Squatina dumeril, Stenotomus chrysops, Syacium papillosum, Symphurus plagiusa, Synagrops bellus, Syngnathus fuscus, Synodus foetens, Synodus poeyi, Tautogolabrus adspersus, Trachinocephalus myops, Trachurus lathami, Trichiurus lepturus, Triglops murrayi, Trinectes maculatus, Ulvaria subbifurcata, Urophycis chuss, Urophycis regia, Urophycis tenuis, Zenopsis conchifera, Zoarces americanus
newf
Alepocephalus bairdii, Amblyraja jenseni, Amblyraja radiata, Ammodytes dubius, Anarhichas denticulatus, Anarhichas lupus, Arctogadus glacialis, Arctozenus risso kroyeri, Argentina silus, Artediellus atlanticus, Artediellus uncinatus, Aspidophoroides olrikii, Bathylagus euryops, Bathyraja richardsoni, Centroscyllium fabricii, Ceramaster granularis, Chiasmodon niger, Chionoecetes opilio, Coryphaenoides guentheri, Cottunculus thomsonii, Crossaster papposus, Eualus gaimardii gaimardii, Eualus pusiolus, Eumicrotremus spinosus, Flabellum macandrewi, Gersemia rubiformis, Glyptocephalus cynoglossus, Gymnelus viridis, Hemitripterus americanus, Hippoglossoides platessoides, Hyas araneus, Icelus bicornis, Illex coindetii, Lebbeus groenlandicus, Leptagonus decagonus, Limanda ferruginea, Liparis atlanticus, Lumpenus fabricii, Lumpenus lampretaeformis, Lycodes lavalaei, Lycodes reticulatus, Malacoraja senta, Mallotus villosus, Merluccius albidus, Micromesistius poutassou, Myoxocephalus aenaeus, Myoxocephalus scorpioides, Myoxocephalus scorpius, Myxine glutinosa, Nemichthys scolopaceus, Notacanthus chemnitzii, Ophiura robusta, Pandalus borealis, Pandalus montagui, Paragorgia arborea, Pasiphaea princeps, Pasiphaea tarda, Pollachius virens, Polyacanthonotus rissoanus, Poromitra capito, Primnoa resedaeformis, Sabinea sarsii, Sabinea septemcarinata, Sebastes norvegicus, Sergestes arcticus, Serrivomer beanii, Spirontocaris phippsii, Stichaeus punctatus, Synaphobranchus kaupii, Triglops murrayi, Trisopterus esmarkii, Urophycis regia
sa
Acanthostracion quadricornis, Achelous spinimanus, Alectis ciliaris, Aluterus schoepfii, Ancylopsetta ommata, Arenaeus cribrarius, Ariopsis felis, Bagre marinus, Bairdiella chrysoura, Brevoortia tyrannus, Calamus leucosteus, Calappa flammea, Callinectes ornatus, Callinectes sapidus, Callinectes similis, Caranx crysos, Carcharhinus acronotus, Carcharhinus limbatus, Caretta caretta, Centropristis philadelphica, Centropristis striata, Chaetodipterus faber, Chilomycterus schoepfii, Chloroscombrus chrysurus, Citharichthys macrops, Citharichthys spilopterus, Cynoscion nothus, Cynoscion regalis, Dasyatis americana, Dasyatis say, Decapterus punctatus, Diplectrum formosum, Dyspanopeus sayi, Echeneis naucrates, Etropus crossotus, Etropus cyclosquamus, Etrumeus teres(sadina), Farfantepenaeus aztecus, Farfantepenaeus duorarum, Gibbesia neglecta, Gymnura micrura, Haemulon aurolineatum, Harengula jaguana, Hepatus epheliticus, Hippocampus erectus, Hypleurochilus geminatus, Lagodon rhomboides, Larimus fasciatus, Leiostomus xanthurus, Litopenaeus setiferus, Lolliguncula brevis, Menippe mercenaria, Menticirrhus americanus, Menticirrhus littoralis, Menticirrhus saxatilis, Micropogonias undulatus, Mobula hypostoma, Mustelus canis, Myliobatis freminvillii, Opisthonema oglinum, Orthopristis chrysoptera, Ovalipes ocellatus, Ovalipes stephensoni, Pagurus pollicaris, Paralichthys albigutta, Paralichthys dentatus, Paralichthys lethostigma, Paralichthys squamilentus, Peprilus paru, Peprilus triacanthus, Persephona mediterranea, Pilumnus sayi, Pomatomus saltatrix, Portunus gibbesii, Prionotus carolinus, Prionotus evolans, Prionotus rubio, Prionotus scitulus, Prionotus tribulus, Rachycentron canadum, Raja eglanteria, Rhinobatos lentiginosus, Rhinoptera bonasus, Rhizoprionodon terraenovae, Sardinella aurita, Scomberomorus cavalla, Scomberomorus maculatus, Scophthalmus aquosus, Selene setapinnis, Selene vomer, Sphoeroides maculatus, Sphyraena guachancho, Sphyrna lewini, Sphyrna tiburo, Squilla empusa, Stellifer lanceolatus, Stephanolepis hispidus, Syacium papillosum, Symphurus plagiusa, Synodus foetens, Trachinotus carolinus, Trichiurus lepturus, Trinectes maculatus, Xiphopenaeus kroyeri
shelf
Alosa pseudoharengus, Alosa sapidissima, Amblyraja radiata, Ammodytes dubius, Anarhichas lupus, Arctozenus risso, Argentina silus, Artediellus atlanticus, Artediellus uncinatus, Aspidophoroides monopterygius, Brosme brosme, Citharichthys arctifrons, Clupea harengus, Cyclopterus lumpus, Enchelyopus cimbrius, Eumicrotremus spinosus, Gadus morhua, Glyptocephalus cynoglossus, Helicolenus dactylopterus, Hemitripterus americanus, Hippoglossoides platessoides, Hippoglossus hippoglossus, Leptagonus decagonus, Leptoclinus maculatus, Leucoraja erinacea, Leucoraja ocellata, Limanda ferruginea, Lophius americanus, Lumpenus lampretaeformis, Lycodes vahlii, Malacoraja senta, Mallotus villosus, Melanogrammus aeglefinus, Merluccius bilinearis, Myoxocephalus octodecemspinosus, Myxine glutinosa, Nezumia bairdii, Peprilus triacanthus, Phycis chesteri, Pollachius virens, Pseudopleuronectes americanus, Reinhardtius hippoglossoides, Scomber scombrus, Squalus acanthias, Triglops murrayi, Urophycis chuss, Urophycis tenuis, Zoarces americanus
wctri
Allosmerus elongatus, Alosa sapidissima, Anoplopoma fimbria, Apristurus brunneus, Argentina sialis, Atheresthes stomias, Bathyagonus pentacanthus, Bathyraja interrupta, Chilara taylori, Citharichthys sordidus, Clupea pallasii pallasii, Cymatogaster aggregata, Engraulis mordax, Enteroctopus dofleini, Eopsetta jordani, Gadus macrocephalus, Genyonemus lineatus, Glyptocephalus zachirus, Gorgonocephalus eucnemis, Hexagrammos decagrammus, Hippasteria phrygiana, Hippoglossoides elassodon, Hippoglossus stenolepis, Hydrolagus colliei, Icelinus filamentosus, Isopsetta isolepis, Loligo opalescens, Lopholithodes foraminatus, Luidia foliolata, Lycodes cortezianus, Lycodes diapterus, Lycodes pacificus, Lyopsetta exilis, Mediaster aequalis, Merluccius productus, Metacarcinus magister, Metridium dianthus, Microgadus proximus, Microstomus pacificus, Oncorhynchus tshawytscha, Ophiodon elongatus, Pandalopsis dispar, Pandalus eous, Pandalus jordani, Pandalus platyceros, Parophrys vetulus, Peprilus simillimus, Pisaster brevispinus, Platichthys stellatus, Pleurobranchaea californica, Pleuronichthys decurrens, Pleuronichthys verticalis, Porichthys notatus, Psettichthys melanostictus, Pycnopodia helianthoides, Raja binoculata, Raja inornata, Raja rhina, Rathbunaster californicus, Rossia pacifica, Scomber japonicus, Sebastes alutus, Sebastes babcocki, Sebastes brevispinis, Sebastes chlorostictus, Sebastes crameri, Sebastes diploproa, Sebastes elongatus, Sebastes entomelas, Sebastes flavidus, Sebastes goodei, Sebastes helvomaculatus, Sebastes jordani, Sebastes paucispinis, Sebastes pinniger, Sebastes proriger, Sebastes reedi, Sebastes ruberrimus, Sebastes saxicola, Sebastes semicinctus, Sebastes wilsoni, Sebastes zacentrus, Sebastolobus alascanus, Squalus acanthias, Strongylocentrotus fragilis, Thaleichthys pacificus, Theragra chalcogramma, Torpedo californica, Trachurus symmetricus, Xeneretmus latifrons, Zalembius rosaceus, Zaniolepis latipinnis Empty data.table (0 rows) of 1 col: reg
naive_msom_scatter()
Figure S1. MSOM richness vs naive richness
MSOM richness and Naive richness are pretty similar. MSOM richness is probably more accurate, or is at least more conservative b/c it has fewer significant trends. But their similarity should help justify using observed presences/ absences in other analyses.
Manuscript paragraph:
MSOM estimates of richness were greater than Naive estimates, but the two methods produced similar temporal dynamics in richness (Figure S1). Henceforth, we report MSOM richness estimates. The greatest long-term average richness was the Gulf of Mexico (142.3), lowest in the Scotian Shelf (46.3). The inter-region average of long-term standard deviations of richness was 5.3; the Aleutian Islands showed the lowest variability (sd = 2.4), and the Gulf of Alaska was the most variable (sd = 8.5).
ur <- comm_master[,unique(reg)]
nr <- length(ur)
naive_msom_mod <- list()
for(r in 1:nr){
td <- comm_master[reg==ur[r]]
naive_msom_mod[[ur[r]]] <- lm(reg_rich~naive_rich, data=td)
}
diff1 <- function(x){
co <- summary(x)$coeff
est <- co["naive_rich","Estimate"]
se <- co["naive_rich","Std. Error"]
# est + sign(est)*1.96*se
pnorm(1-est, sd=se, lower.tail=FALSE)
}
rbindlist(lapply(naive_msom_mod, mod_smry, pred_name="naive_rich"), idcol='reg')
## reg Class mod_call predictor estimate p.value Marginal
## 1: ai lm reg_rich ~ naive_rich naive_rich 0.76768 0e+00 0.9094
## 2: ebs lm reg_rich ~ naive_rich naive_rich 0.85367 0e+00 0.9799
## 3: gmex lm reg_rich ~ naive_rich naive_rich 0.70116 0e+00 0.8240
## 4: goa lm reg_rich ~ naive_rich naive_rich 0.76002 2e-04 0.7283
## 5: neus lm reg_rich ~ naive_rich naive_rich 1.04916 0e+00 0.8834
## 6: newf lm reg_rich ~ naive_rich naive_rich 1.06926 0e+00 0.9917
## 7: sa lm reg_rich ~ naive_rich naive_rich 0.83276 1e-03 0.3797
## 8: shelf lm reg_rich ~ naive_rich naive_rich 0.74526 0e+00 0.9006
## 9: wctri lm reg_rich ~ naive_rich naive_rich 0.97977 0e+00 0.9914
## Conditional AIC
## 1: NA 30.333
## 2: NA 78.322
## 3: NA 73.003
## 4: NA 80.570
## 5: NA 161.988
## 6: NA 20.860
## 7: NA 133.744
## 8: NA 123.260
## 9: NA 25.122
sapply(naive_msom_mod, diff1)
## ai ebs gmex goa neus newf
## 1.2144e-03 5.9555e-11 1.7715e-04 4.3195e-02 7.6001e-01 9.9597e-01
## sa shelf wctri
## 2.2556e-01 6.5869e-11 2.6534e-01
All regions except Northeast U.S. and Newfoundland have estimates that are less than 1, although West Coast U.S. is not signifiantly less than 1. Therefore, for two-thirds of the regions, the MSOM tended to estimate that more undetected species were present when naive estimates were low.
sumry_counts <- data_all[reg!="wcann",
j = {
yi <- table(diff(sort(unique(year))))
yi_form <- paste0(names(yi), paste("(",yi,")",sep=""))
data.table(
"N" = trawlData::lu(year),
"Years" = paste(min(year),max(year),sep=" - "),
"Frequency" = paste(yi_form,collapse=", "),
"Sites" = trawlData::lu(stratum),
"Tows" = paste0(
.SD[,trawlData::lu(haulid),by=c("stratum","year")][,max(V1)],
paste0("(",.SD[,trawlData::lu(haulid),by=c("stratum","year")][,round(mean(V1),1)],")")
),
"Species" = trawlData::lu(spp)
)
},
by=c('reg')
]
setnames(sumry_counts, 'reg', "Region")
sumry_counts[,Region:=factor(Region, levels=c("ebs","ai","goa","wctri","gmex","sa","neus","shelf","newf"), labels=c("E. Bering Sea","Aleutian Islands","Gulf of Alaska","West Coast US","Gulf of Mexico","Southeast US", "Northeast US", "Scotian Shelf","Newfoundland"))]
setkey(sumry_counts, Region)
sumry_counts[,Org:=c("AFSC","AFSC","AFSC","NMFS","GSMFC","SCDNR","NEFSC","DFO","DFO")]
stargazer(sumry_counts, summary=FALSE, rownames=FALSE, column.sep.width="2pt", digits.extra=0, digits=NA)
% Table created by stargazer v.5.2 by Marek Hlavac, Harvard University. E-mail: hlavac at fas.harvard.edu % Date and time: Wed, Jun 28, 2017 - 15:54:04
# \begin{table}[!htbp] \centering
# % \caption{}
# % \label{}
# \begin{tabular}{@{\extracolsep{2pt}} cccccccc}
# \\[-1.8ex]\hline
# \hline \\[-1.8ex]
# Region & N & Years & Frequency & Sites & Tows & Species & Org \\
# \hline \\[-1.8ex]
# E. Bering Sea & $31$ & 1984 - 2014 & 1(30) & $138$ & 4(1.6) & $110$ & $\text{AFSC}^{*}$ \\
# Aleutian Islands & $12$ & 1983 - 2014 & 2(5), 3(4), 4(1), 5(1) & $82$ & 12(2.8) & $54$ & $\text{AFSC}^{*}$ \\
# Gulf of Alaska & $13$ & 1984 - 2013 & 2(7), 3(5) & $89$ & 15(4.4) & $98$ & $\text{AFSC}^{*}$ \\
# West Coast US & $10$ & 1977 - 2004 & 3(9) & $84$ & 91(4.7) & $92$ & $\text{NMFS}^{\dagger}$ \\
# Gulf of Mexico & $17$ & 1984 - 2000 & 1(16) & $39$ & 39(5.7) & $144$ & $\text{GSMFC}^{\ddagger}$ \\
# Southeast US & $25$ & 1990 - 2014 & 1(24) & $24$ & 13(3.9) & $104$ & $\text{SCDNR}^{\S}$ \\
# Northeast US & $32$ & 1982 - 2013 & 1(31) & $100$ & 10(3) & $141$ & $\text{NEFSC}^{\P}$ \\
# Scotian Shelf & $41$ & 1970 - 2010 & 1(40) & $48$ & 11(2.5) & $48$ &$\text{DFO}^{\nabla}$ \\
# Newfoundland & $16$ & 1996 - 2011 & 1(15) & $191$ & 9(2.2) & $72$ & $\text{DFO}^{\nabla}$ \\
# \hline \\[-1.8ex]
# \end{tabular}
# \end{table}
# knitr::kable(sumry_counts, caption="Years = Total number of years sampled, Year Range = the minimum and maximum of years sampled, Year Interval = the number of years elapsed between samples (and the frequency of this interval in parentheses), Sites = the total number of sites in the region, Max. Tows = maximum number of tows per site per year, Average Tows = the average number of tows per site per year, Total Species = the total number of species observed in the region across all sites and years.")
| reg | estimate | BH | pvalue |
|---|---|---|---|
| ai | -0.09091 | 0.73170 | 0.73170 |
| ebs | 0.50108 | 0.00018 | 0.00008 |
| gmex | 0.85294 | 0.00001 | 0.00000 |
| goa | 0.20513 | 0.40514 | 0.36012 |
| neus | 0.33681 | 0.01032 | 0.00803 |
| newf | 0.78658 | 0.00011 | 0.00004 |
| sa | -0.56854 | 0.00028 | 0.00016 |
| shelf | 0.70488 | 0.00000 | 0.00000 |
| wctri | 0.76409 | 0.00456 | 0.00304 |
# par(cex=1, mar=c(3,1,1,0), ps=8, mgp=c(1,0.25,0), tcl=-0.15, lheight=0.7)
categ_barplot()
# png("~/Desktop/bp.png", width=3.5, height=3.5, res=200, units='in'); categ_barplot(); dev.off()
Figure S2. Number of species beloning to the categories of both, neither, colonizer, leaver in each region
| ai | ebs | gmex | goa | neus | newf | sa | shelf | wctri | |
|---|---|---|---|---|---|---|---|---|---|
| neither | 39 | 64 | 105 | 63 | 56 | 49 | 69 | 27 | 64 |
| both | 6 | 44 | 31 | 17 | 83 | 12 | 35 | 20 | 15 |
| colonizer | 8 | 2 | 8 | 18 | 1 | 10 | 0 | 1 | 12 |
| leaver | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 1 |
It’s the same pattern, whichever way you split it. However, AI is the only region that had more colonizers than both species. An interesting way to think about some of this is that the average sd in richness was 5.29901, so when the number of colonizer or leaver species exceed’s that region’s sd, the impact of those categories, which I consider to be dubious, might start being relevant (though it’s not necessarily problematic, nor is this even close to an actual test for the significance of those categories to the trend). EBS and Shelf had significant positive trends in richness and very low numbers in the colonizer category. WCTRI and NEWF had similar numbers in the both and colonizer category.
| Region | Richness Change | colonizer |
|---|---|---|
| E. Bering Sea | 12.4666 | 2 |
| Aleutian Islands | 4.2513 | 8 |
| Gulf of Alaska | 9.5026 | 1 |
| West Coast US | 7.3750 | 8 |
| Gulf of Mexico | -5.0019 | 0 |
| Southeast US | 18.8206 | 12 |
| Northeast US | 12.7363 | 10 |
| Scotian Shelf | 7.5842 | 1 |
| Newfoundland | 17.5343 | 18 |
stargazer(att_categ_print, summary=FALSE, rownames=FALSE, column.sep.width="2pt", digits=1, digits.extra=0) # print table
% Table created by stargazer v.5.2 by Marek Hlavac, Harvard University. E-mail: hlavac at fas.harvard.edu % Date and time: Wed, Jun 28, 2017 - 15:54:05
# \begin{table}[!htbp] \centering
# %\caption{}
# %\label{}
# \begin{tabular}{@{\extracolsep{2pt}} ccc}
# \\[-1.8ex]\hline
# \hline \\[-1.8ex]
# Region & Richness Change & Colonizing Species \\
# \hline \\[-1.8ex]
# E. Bering Sea & $12.5$ & $2$ \\
# Aleutian Islands & $4.3$ & $8$ \\
# Gulf of Alaska & $17.5$ & $18$ \\
# West Coast US & $18.8$ & $12$ \\
# Gulf of Mexico & $7.4$ & $8$ \\
# Southeast US & $$-$5.0$ & $0$ \\
# Northeast US & $9.5$ & $1$ \\
# Scotian Shelf & $7.6$ & $1$ \\
# Newfoundland & $12.7$ & $10$ \\
# \hline \\[-1.8ex]
# \end{tabular}
# \end{table}
ceEventRange("mean_size")
Figure S7. Number of colonizations and extinctions as a function of range size and range density.
Yup, this is definitely a thing. Long-term average range size and range density predict how many colonizations and extinctions a species is likely to have. This will lead nicely into examing how range changes prior to an extinction or after a colonization.
par(mfrow=c(3,3), mgp=c(0.85,0.2,0), ps=8, cex=1, mar=c(1.75,1.5,0.75,0.5), oma=c(0.5,0.5,0.1,0.1), tcl=-0.15)
ur <- names(pretty_reg)
for(r in 1:length(ur)){
t_reg <- ur[r]
t_dat <- rangeTimeDT[reg==t_reg]
t_dat[,spp:=paste0(spp,event)]
t_dat[type=="pre_ext", c("spp","time"):=list(spp=paste0(spp,type),time=-time)]
setkey(t_dat, spp, type, time)
scatterLine(t_dat, x="time", y="size", lineBy="spp", colBy=adjustcolor(pretty_col[t_reg], 1), lwdBy=0.5, type='p', xlab="", ylab="")
if(t_reg=="ai"){
t_dat[time <=-10,lines(time, size, col='black')]
}
t_dat[,points(mean(time), mean(size), pch=19, col='black', cex=1.2),by=c('time')]
t_dat[,points(mean(time), mean(size), bg=pretty_col[t_reg], pch=21, col='white'),by=c('time')]
mtext(pretty_reg[t_reg], line=0, side=3, adj=0.1, font=2)
}
mtext("Time to Event", side=1, line=-0.5, outer=TRUE)
mtext("Range Size", side=2, line=-0.5, outer=TRUE)
col_ext_ts()
Figure NotIncluded. Number of colonizations (blue) and extinctions (red) over time in each region.
In most regions the differences in colonization and extinction numbers are similar over time. The most obvious exceptions are for the 3 regions that showed large initial spikes in richness; the GOA, GMEX, and AI regions initially have much larger numbers of colonizers than leavers, but this number shrinks rapidly until the two rates are ~equal.
For the regions with significant positive slopes, there is no visually obvious increase in colonizations relative to extinctions over time. Because the colonization and extinction numbers tend to track each other over the long-term, it it would be difficult to attribute the long-term changes in richness to a change in just colonization or extinction rates.
Manuscript paragraph:
A time series of richness can be decomposed into the colonizations and extinctions of individual species over time. We categorized species according to the following colonization extinction patterns: present in all years = neither (536 species), colonized and went extinct = both (263 species), initially absent but present every year after its colonization = colonizer (61 species), initially present but absent every year after its extinction = leaver (4 species). Most regions had the same overall ranking (neither > both > colonizer > leaver), except in the Northeast US where both was the most common and neither second, and in the Aleutian Islands where colonizer was the second most common and both third (Figure S2). In general, changes in richness were not due to permanent departures or introductions of species to the region. Furthermore, colonization and extinction rates did not become more dissimilar over time for any region (Figure NotIncluded). Colonizations were initially greater than extinctions in Aleutian Islands, Gulf of Alaska, and Gulf of Mexico, but this difference disappeared in the latter portion of these time series, as evidenced by these regions’ initially rapid increase in richness that later plateaued. The other regions did not show strong trends in the difference between colonizations and extinctions over time, making it difficult to attribute the long-term trends in richness to changes in just colonization or just extinction rates.
rangeSize_absenceTime("rangeDensity")
Figure 4b. Geographic range density vs years until extinction (A) and years after colonization (B). For visualization purposes, range density is averaged across species for each unique value on each axis, and a linear model fit through this average. Statistics in main text use unaggregated data. The horizontal axes were formulated as time until (since) the nearest upcoming (previous) absence. Because range density must be zero when either horizontal axis has a value of zero, points at (0,0) were excluded from figures and analyses.
fp <- spp_master[present==1,.SD[which.min((year))] ,by=c('reg','spp')]
fp2 <- fp[,list(new_col=sum(col)),keyby=c('reg','year')]
ry_skele <- unique(spp_master[,list(reg,year)])
fp3 <- merge(ry_skele, fp2, all=TRUE, by=c('reg','year'))
fp3[is.na(new_col), new_col:=0]
par(mfrow=c(3,3))
fp3[,plot(year, new_col, type='o', main=reg[1]),by=c('reg')]
Empty data.table (0 rows) of 1 col: reg
Reveals similar information as the number of colonizations per year, except it isolates the number of colonizations by species that have never been seen in the region before. This is the rate at which new species are introduced to the region. An important caveat here is that the trends here are probably biased downward because we onlly included species that were observed in at leat 10 year-sites in order to avoid the inclusion of ultra-rare species. These exclusions reduce the odds that species that colonized a limited number of sites towards the end of the time series met the threshold of 10, especially for regions with relatively few strata.
E.g., say that a species colonized Southeast U.S., which has 24 sites, in the last year of sampling — for this species to be included, it would have had to occupy 42% of the region (10 sites) in its first year present. Across all regions, the average colonizing species occupies only 5% of strata in its first year, and for the 3 species that colonized the Southeast U.S., they colonized only 4.2% of sites, a full order of magnitude lower than what was needed. Depending on the expansion rate, a species would need to first colonize several years before the end of the time series to be included. It’s also worth noting that no region has a new colonizing species towards the end of the time series — this makes a similar point.
Because of the way we eliminated species from our data set, our analysis was extremely conservative with regard to detecting changes in richness, particularly if those changes were going to be the result of changes in the rate of novel introductions (right now I’m a little less clear on what sources of species declines might be affected, other than to say that generally fewer species would have been included due to only being present at the beginning [affects declines more] or end [affects increases more] of the time series). Regions with short time series and few sites were especially prone to having muted changes in richness.
rangeSizeDens()
Figure S6. Range density versus range size. In panel A each point is a species-region-year combination. In panel B, each point is a region-year. Range size is the proportion of sites occupied, range density the tows in occupied sites. The community metrics in B is calculated by take each species’ long-term average from A, then taking the average across all species present in the community in a given year. Fitted lines in A are from a loess fit.
First posterity, I’ll show a figure relating range size and range density. Then I’ll present the figure of how size can ‘predict’ richness. Finally, I’ll explore models more formally expressing the relationship between richness and size.
These plots show that there is a strong relationship between range size and range density. Interestingly, in Figure S6B, the cross-region relationship is negative (if each color had 1 pt), whereas the within-region relationship is positive.
The positive relationship between size and density is not surprising. My interpretation of density is ~population size. Often population size is correlated with range size. I think this is a standard result, but I need to double-check.
##Figure. Range size vs absence time as individual regressions
rangeTimeDT[,nTime:=length(time),by=c("reg","type","event","spp")]
o <- rangeTimeDT[nTime>=3,j={
getEPR(lm(size~time))
} ,by=c("reg","type","event","nTime","spp")
]
rangeTime_signif <- o[,list(propSignificant=(sum(Pr<0.05)/sum(!is.na(Pr))),n=sum(!is.na(Pr))),by=c("reg","type")]
# histograms of estimates, p-value, and r-squared split by pre-extinction and post-colonization
par(mfrow=c(4,2), mar=c(2,2,0.5,0.5), cex=1, ps=10, mgp=c(1,0.1,0), tcl=-0.1)
o[,j={hist(Estimate, main=type[1]);NULL},by='type']
Empty data.table (0 rows) of 1 col: type
o[,j={hist(Pr, main=type[1]);NULL},by='type']
Empty data.table (0 rows) of 1 col: type
o[,j={hist(Rsquared, main=type[1]);NULL},by='type']
Empty data.table (0 rows) of 1 col: type
# p value vs number of years in series
setorder(o, nTime)
# o[,j={plot(nTime,Pr); lines(spline(Pr~nTime),col='red',lwd=2)}]
# abline(h=0.05, lwd=1.5)
# abline(h=0.05, col='white', lwd=0.5)
o[,mean(Pr),by='nTime'][,plot(nTime, V1)];
NULL
abline(h=0.05, col='blue', lty=2)
# slope vs number of years in series
# o[,j={plot(nTime,Estimate); lines(spline(Estimate~nTime),col='red',lwd=2)}]
# abline(h=0, lwd=1.5)
# abline(h=0, col='white', lwd=0.5)
o[,mean(Rsquared),by='nTime'][,plot(nTime, V1)];
NULL
Exploration Figure histograms of separate regressions of size ~ time; this is for each run-up to an extinction and reach follow-up to a colonization. Trying to understand how regularly the pattern might be observed. Hard to answer because adjusting the restriction on number of events in the run-up or follow-up (nTime) greatly affects the proportion that are significant.
| reg | type | propSignificant | n |
|---|---|---|---|
| ai | post_col | 0.35714 | 14 |
| ai | pre_ext | 0.00000 | 2 |
| ebs | pre_ext | 0.27027 | 37 |
| ebs | post_col | 0.35294 | 51 |
| gmex | post_col | 0.16216 | 37 |
| gmex | pre_ext | 0.08333 | 12 |
| goa | post_col | 0.25714 | 35 |
| goa | pre_ext | 0.25000 | 4 |
| neus | post_col | 0.17647 | 102 |
| neus | pre_ext | 0.05714 | 70 |
| newf | post_col | 0.21053 | 19 |
| newf | pre_ext | 0.16667 | 6 |
| sa | pre_ext | 0.27027 | 37 |
| sa | post_col | 0.08571 | 35 |
| shelf | post_col | 0.30769 | 26 |
| shelf | pre_ext | 0.06667 | 15 |
| wctri | post_col | 0.21739 | 23 |
| wctri | pre_ext | 0.50000 | 2 |
| type | proportionSignificant |
|---|---|
| post_col | 0.23635 |
| pre_ext | 0.18493 |
The mixed effect models show a strong tendency for range size to be smaller near an absence. I figured that this relationship wouldn’t be apparent for all cases, which ended up being the result. However, the significance of these relationships depends on the number of years for each event stretch. So sample size matters. On one hand, this could hold implications for the rate of decline, so excluding short stretches might also exclude more cases with abrupt declines/ increases (although, these may well be preceded by long durations of stability, so not all abrupt declines/ increases would be excluded necessarily). However, the slope estimate doesnt change a whole lot across sample size: from 0 to ~10 years, increasing sample size (duration) also increases the slope (going from slightly negative to pretty consistently positive). Then it levels out, and doesn’t continue increasing. So it is not the case that longer stretches are long because they have more gradual slopes, necessarily. Actually, in the range of data for which there does appear to be a trend (0-10 samples), the slope becomes larger which is the opposite of “longer stretches result from more gradual processes” notion (thinking in terms of extinction).
Ultimately, I think the mixed effects model is far more appropriate so long as conclusions are cauched in general patterns, not in being able to detect the pattern for any single species. While it could have been interesting to present the many-regression results too, I think it is a separate analysis to consider how often this rule would apply to individual species. Such an analysis would benefit from understanding when the rule does and does not apply, not just the apparent frequency of relevance.
blah <- spp_master[present==1, j={
list(propSlope=summary(lm(range_size_mu~year))$coeff[2,1], mu_propStrat=mean(range_size_mu), ce_categ=ce_categ[1])
} ,by=c("reg","spp")]
par(mfrow=c(3,3))
ureg <- blah[,unique(reg)]
nreg <- length(ureg)
for(r in 1:nreg){
blah[reg==ureg[r],j={boxplot(propSlope~ce_categ, main=reg[1], outline=FALSE);abline(h=0);NULL}]
}
I wanted to determine if species in the region are becoming more widespread over time. In other words, when the species is present, what fraction of the region does it occupy? Does that fraction change over time? This slope is the y-axis for the boxplots. It indicates the long-term slope of the proportion of sites occupied (excluding any year when that proportion was 0). The x-axis represents the categories of colonization-extinction patterns. Species in the “both” category experienced both colonization and extinction events at some point in the time series.
Other figures and analyses suggest that species in the “both” category play an important role in increases in species richness. In other words, the coming-and-going of these organisms is not neutral in the long term. This conjecture is supported by 1) the magnitude of increase for regions with increasing richness cannot be explained by “colonizing” species alone; 2) geographically constrained species have more colonizations and extinctions that widespread species.
Looking at the different box categories, if historically rare species (i.e., those that were in the both category and had most of the colonizations/ extinctions) were appearing more consistently in the time series because they were becoming more geogrpahically widespread (and therefore less likely to go extinct), then we’d expect to see that the “both” category to have more-positive slopes in their annual range sizes.
That’s kinda? what I see. It’s definitely not a very strong relationship. One issue with this analysis is that you wouldn’t expect any given species in the “both” category to become more prevalent. In fact, it should only be a small fraction of those species that are becomign more widespread and thus contributing to long-term increases in species richness.
localR <- data_all[,list(lR=length(unique(spp))),by=c("reg","stratum","year")]
bothR <- merge(localR, comm_master[,list(reg,year,naive_rich,reg_rich)],by=c("reg","year"))
par(mfrow=c(3,3), cex=1, ps=8, mar=c(2,2,2,2), mgp=c(1,0.2,0), tcl=-0.2)
ureg <- bothR[,unique(reg)]
nreg <- length(ureg)
for(r in 1:nreg){
comm_master[reg==ureg[r],j={
plot(year, get(lR), type='l', ylab=lR)
mtext(reg[1], side=3, line=0, adj=0.1, font=2)
par(new=TRUE)
plot(year, get(rR), type='l', col='red', xaxt='n', yaxt='n', xlab='',ylab='')
axis(side=4, col='red')
mtext(rR, side=4, line=1, col='red')
NULL
}]
}
mtext("Observed Regional (red) and Mean Local (black) Richness", side=3, line=-0.75, outer=TRUE, font=2)
In general, local richness and regional richness are similar. There are differences, possibly in some cases the regional slope might be significant whereas local would not be (I haven’t checked, though). However, the trends aren’t in opposite directions at the two scales.
eval(figure_setup())
par(mfrow=c(3,1), mar=c(1.75,1.5,0.25,0.25),mgp=c(0.85,0.1,0), tcl=-0.1, cex=1, ps=8)
ureg <- comm_master[,unique(reg)]
nreg <- length(ureg)
ptCol <- bquote(adjustcolor(pretty_col[reg],0.5))
scatterLine(Data=comm_master, x="naive_rich", y="beta_div_obs", lineBy="reg", col=ptCol, pch=20)
scatterLine(Data=comm_master, x="local_rich_obs", y="beta_div_obs", lineBy="reg", col=ptCol, pch=20)
comm_master[,legend("topright",ncol=2,legend=pretty_reg[una(reg)],text.col=pretty_col[una(reg)], inset=c(-0.01, -0.03), bty='n', x.intersp=0.15, y.intersp=0.65)]
rect text
1: 32.124 30.318,30.318,30.318,30.318,30.318,45.775, 2: 0.081876 0.39137,0.37885,0.36632,0.35380,0.34128,0.39137, 3: 27.532 30.318,30.318,30.318,30.318,30.318,45.775, 4: 0.40726 0.39137,0.37885,0.36632,0.35380,0.34128,0.39137,
scatterLine(Data=comm_master, x="naive_rich", y="local_rich_obs", lineBy="reg", col=ptCol, pch=20)
eval(figure_setup())
par(mfrow=c(3,1), mar=c(1.75,1.5,0.25,0.25),mgp=c(0.85,0.1,0), tcl=-0.1, cex=1, ps=8)
ureg <- comm_master[,unique(reg)]
nreg <- length(ureg)
ptCol <- bquote(adjustcolor(pretty_col[reg],0.5))
scatterLine(Data=comm_master, x="reg_rich", y="beta_div_mu", lineBy="reg", col=ptCol, pch=20)
scatterLine(Data=comm_master, x="local_rich", y="beta_div_mu", lineBy="reg", col=ptCol, pch=20)
comm_master[,legend("topright",ncol=2,legend=pretty_reg[una(reg)],text.col=pretty_col[una(reg)], inset=c(-0.01, -0.03), bty='n', x.intersp=0.15, y.intersp=0.65)]
rect text
1: 50.059 48.689,48.689,48.689,48.689,48.689,72.775, 2: 0.11305 0.36096,0.34367,0.32638,0.30909,0.29180,0.36096, 3: 44.348 48.689,48.689,48.689,48.689,48.689,72.775, 4: 0.38291 0.36096,0.34367,0.32638,0.30909,0.29180,0.36096,
scatterLine(Data=comm_master, x="reg_rich", y="local_rich", lineBy="reg", col=ptCol, pch=20)
noNeither <- spp_master[,j={
td <- .SD[ce_categ!="neither" & present==1]
td1.0 <- td[,list(propStrata_noNeither_ltAvg=mean(range_size_mu)),keyby=c("reg","spp")]
setkey(td, reg, spp)
td1.1 <- td[td1.0]
td1.2 <- td1.1[,list(propStrata_noNeither_avg_ltAvg=mean(propStrata_noNeither_ltAvg)),by=c("reg","year")]
td2 <- td[,list(propStrata_noNeither_avg=mean(range_size_mu)),by=c("reg","year")]
td_out <- merge(td1.2, td2, by=c('reg','year'))
}]
cmNN <- merge(comm_master, noNeither)
par(mfrow=c(2,2))
ureg <- comm_master[,unique(reg)]
nreg <- length(ureg)
ptCol <- bquote(adjustcolor(pretty_col[reg],0.5))
scatterLine(Data=cmNN, x="range_size_mu_avg_ltAvg", y="reg_rich", lineBy="reg", col=ptCol, pch=20)
scatterLine(Data=cmNN, x="range_size_mu_avg", y="reg_rich", lineBy="reg", col=ptCol, pch=20)
scatterLine(Data=cmNN, x="propStrata_noNeither_avg_ltAvg", y="reg_rich", lineBy="reg", col=ptCol, pch=20)
scatterLine(Data=cmNN, x="propStrata_noNeither_avg", y="reg_rich", lineBy="reg", col=ptCol, pch=20)
std_range_t1 <- spp_master[ce_categ!='neither' & ce_categ!="colonizer", j={
std_range <- .SD[,list(year, ext, std_range=c((get(rSMet_base))), ext_dist=ext_dist),by=c("spp")]
setkey(std_range, spp, year)
std_range_t1 <- std_range[present==1,list(year=year[-1], ext=ext[-1], std_range=std_range[-1], std_range_t1=head(std_range, -1), ext_dist=ext_dist[-1]),by=c("spp")]
std_range_t1[,nspp:=length(unique(spp)),by=c('year')]
std_range_t1
},by="reg"]
# print(std_range_t1[,{print(summary(glm(ext~std_range, family='binomial')));NULL},by='reg'])
# print(std_range_t1[,{print(summary(glmer(ext~std_range+(std_range|reg), family='binomial')));NULL}])
# Set up -- region names and empty named lists
ur <- std_range_t1[,unique(reg)]
rangeExtProb1 <- structure(vector("list",length(ur)), .Names=ur)
rangeExtProb2 <- structure(vector("list",length(ur)), .Names=ur)
rangeExtProb3 <- structure(vector("list",length(ur)), .Names=ur)
# Fit 3 models for each region
for(r in 1:length(ur)){
t_reg <- ur[r]
t_dat <- std_range_t1[reg==t_reg]
rangeExtProb1[[t_reg]] <- (glm(ext~std_range, family='binomial', data=t_dat))
rangeExtProb2[[t_reg]] <- summary(glm(ext~std_range+nspp, family='binomial', data=t_dat))
rangeExtProb3[[t_reg]] <- summary(glm(ext~std_range*nspp, family='binomial', data=t_dat))
}
rangeExtProb1_smry <- (smry_modList(rangeExtProb1, pred_name="std_range"))[,reg:=names(rangeExtProb1)]
rangeExtProb1_smry[,BH:=p.adjust(p.value, method='BH')]
par(mfrow=c(3,3), mar=c(2.5,2.5,1.5,0.5), ps=10, mgp=c(0.75,0.15,0), tcl=-0.15)
for(r in 1:length(ur)){
std_range_t1[reg==ur[r],plot(std_range, ext)]
new_r <- std_range_t1[reg==ur[r],sort(std_range)]
lines(new_r, predict(rangeExtProb1[[ur[r]]], newdata=data.frame(std_range=new_r), type='response'))
mtext(pretty_reg[ur[r]], side=3, line=0, adj=0.1, font=2)
}
make_timeRange <- function(densName=c("propTow_occ"), sizeName=c("range_size_samp","range_size_mu","propStrata")){
rangeTimeDT <- spp_master[stretch_type=="pre_ext" & !is.na(stretch_type) & propStrata!=0]
rangeTimeDT <- rangeTimeDT[,list(
reg=reg,
event=as.character(event_year),
spp=spp,
type=as.character(stretch_type),
time=ext_dist,
size=get(sizeName),
density=get(densName)
)]
return(rangeTimeDT)
}
timeRangeDT <- make_timeRange()
ur <- timeRangeDT[,unique(reg)]
sTime_reg_mods_Inv <- structure(vector("list",length(ur)), .Names=ur)
for(r in 1:length(ur)){
t_reg <- ur[r]
t_dat <- timeRangeDT[reg==t_reg]
sTime_reg_mods_Inv[[t_reg]] <- lme4::lmer(time ~ size + (size|spp), data=t_dat)
}
sTime_reg_smry_Inv <- smry_modList2(sTime_reg_mods_Inv)
kable(sTime_reg_smry_Inv)
| reg | mod_call | MargR2 | CondR2 | AIC | pval_size | BH_size | randGrp | Int | size |
|---|---|---|---|---|---|---|---|---|---|
| ai | time~size+(size|spp) | 0.10570 | 0.2884 | 173.60 | 0.09039 | 0.203 | spp | 6.7767 | 13.7420 |
| ebs | time~size+(size|spp) | 0.19780 | 0.4588 | 1773.99 | 0.00000 | 0.000 | spp | 2.5463 | 29.0803 |
| gmex | time~size+(size|spp) | 0.00764 | 0.3658 | 641.43 | 0.47281 | 0.473 | spp | 2.8071 | 10.4723 |
| goa | time~size+(size|spp) | NA | NA | 179.91 | 0.19779 | 0.297 | spp | 3.5593 | 53.2009 |
| neus | time~size+(size|spp) | 0.00505 | 0.5083 | 4075.28 | 0.41940 | 0.472 | spp | 3.1922 | 15.3583 |
| newf | time~size+(size|spp) | 0.03149 | 0.4440 | 351.84 | 0.16245 | 0.292 | spp | 3.5170 | 8.4463 |
| sa | time~size+(size|spp) | 0.23210 | 0.5018 | 1969.26 | 0.00000 | 0.000 | spp | 3.0907 | 23.9555 |
| shelf | time~size+(size|spp) | 0.08943 | 0.3414 | 648.72 | 0.02742 | 0.082 | spp | 2.2084 | 33.8794 |
| wctri | time~size+(size|spp) | 0.12500 | 0.7858 | 153.90 | 0.26425 | 0.340 | spp | 3.6162 | 49.8260 |
par(mfrow=c(3,3), mar=c(1.85,1.85, 0.5,0.5), mgp=c(0.65, 0.15, 0), oma=c(0.5, 0.5, 0.25, 0.1), tcl=-0.15, ps=8, cex=1)
ur <- names(pretty_reg)[names(pretty_reg)%in%comm_master[,unique(reg)]]
nr <- length(ur)
for(r in 1:nr){
comm_master[reg==ur[r],j={
plot(beta_div_obs, propStrata_avg_ltAvg, xlab='', ylab='')
mod <- lm(propStrata_avg_ltAvg~beta_div_obs)
abline(mod)
}]
mtext(pretty_reg[ur[r]], side=3, adj=0.1, font=2)
}
eval(figure_setup())
eval(neitherPres_subColor())
boxRange_colRich(range_type="range_size_samp")
dat <- merge(transExpansionStats, comm_master[,coef(lm(reg_rich~year))[2], by='reg'], by='reg')
dat[,summary(lm(V1~year))]
Call: lm(formula = V1 ~ year)
Residuals: Min 1Q Median 3Q Max -0.3989 -0.1088 0.0060 0.0928 0.3423
Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.170 0.111 1.53 0.169
year 106.791 38.768 2.75 0.028 * — Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘.’ 0.1 ‘’ 1
Residual standard error: 0.238 on 7 degrees of freedom Multiple R-squared: 0.52, Adjusted R-squared: 0.452 F-statistic: 7.59 on 1 and 7 DF, p-value: 0.0283
dat[,plot(year, V1, col=pretty_col[reg], pch=19, xlab="Rate of Range Expansion by Transient Species", ylab="Rate of Change in Species Richness")]
NULL
The significance of this trend depends on Southeast US, though. I’m not sure I want to get into defending a plot based on an outlier, and the fact that I didn’t account for uncertainty in the x-axis.
Sys.time()
## [1] "2017-06-28 15:54:16 EDT"
sessionInfo()
## R version 3.3.0 (2016-05-03)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X 10.11.6 (El Capitan)
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] texreg_1.36.4 stargazer_5.2 kfigr_1.2
## [4] xtable_1.8-2 rmarkdown_1.0 knitr_1.14
## [7] viridis_0.4.0 viridisLite_0.2.0 piecewiseSEM_1.1.3
## [10] spdep_0.6-5 spatstat_1.45-2 rpart_4.1-10
## [13] nlme_3.1-127 raster_2.5-2 sp_1.2-3
## [16] maps_3.1.0 multcomp_1.4-5 TH.data_1.0-7
## [19] MASS_7.3-45 survival_2.39-2 mvtnorm_1.0-5
## [22] car_2.1-2 lme4_1.1-12 Matrix_1.2-6
## [25] rbLib_0.1.1 trawlDiversity_0.5.0 trawlData_0.4.0
## [28] data.table_1.10.4
##
## loaded via a namespace (and not attached):
## [1] bit64_0.9-5 splines_3.3.0 gtools_3.5.0
## [4] highr_0.6 stats4_3.3.0 yaml_2.1.13
## [7] LearnBayes_2.15 pbivnorm_0.6.0 lattice_0.20-33
## [10] quantreg_5.21 quadprog_1.5-5 digest_0.6.10
## [13] polyclip_1.5-6 ggm_2.3 minqa_1.2.4
## [16] colorspace_1.2-6 sandwich_2.3-4 htmltools_0.3.5
## [19] plyr_1.8.4 SparseM_1.7 gmodels_2.16.2
## [22] scales_0.4.0 gdata_2.17.0 tensor_1.5
## [25] MatrixModels_0.4-1 mgcv_1.8-12 ggplot2_2.1.0
## [28] nnet_7.3-12 pbkrtest_0.4-6 mnormt_1.5-4
## [31] magrittr_1.5 deldir_0.1-12 evaluate_0.9
## [34] tools_3.3.0 formatR_1.4 stringr_1.2.0
## [37] munsell_0.4.3 grid_3.3.0 nloptr_1.0.4
## [40] goftest_1.0-3 igraph_1.0.1 lavaan_0.5-20
## [43] boot_1.3-18 gtable_0.2.0 codetools_0.2-14
## [46] abind_1.4-3 gridExtra_2.2.1 zoo_1.7-13
## [49] bit_1.1-12 stringi_1.1.5 parallel_3.3.0
## [52] Rcpp_0.12.10 coda_0.18-1